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fluid dynamics equations 

Y. Saad * D. Semeraro t 


Abstract 

In this paper we present an application of matrix exponentiation via Krylov sub- 
space projection, to the solution of fluid dynamics problems. The main idea is to ap- 
proximate the operation exp(A)v by means of a projection-like process onto a Krylov 
subspace. This results in a computation of an exponential matrix vector product simi- 
lar to the one above but of a much smaller size. Time integration schemes can then be 
devised to exploit this basic computational kernel. The motivation of this approach is 
to provide time-integration schemes that are essentially of an explicit nature but which 
have good stability properties. 
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1 Introduction 

This paper presents some results of the application of a new integration technique to the 
Navier-Stokes equations. The originality of the new method is that it attempts to exploit 
the local exponential behavior of the solution with respect to time. Roughly speaking, 
evaluating the solution of a nonlinear system of equation w f = F^w) at time t + A t involves 
the computation of a vector of the form exp(AtA)w(t), where w(t) is the approximation at 
the current time and A is the linearization of the nonlinear mapping F at w(t). This type 
of approximation is exact for linear problems. Although it may seem difficult to exploit 
this local bahavior in practice because it involves the computation of the exponential of a 
matrix times a vector, we will show that using a Krylov subspace approach, a number of 
efficient integration schemes can be derived based on it. Moreover the scheme is of an explicit 
nature in that the only operations that it requires with the Jacobian A are matrix by vector 
products. 

Many of the integration schemes for solving time dependent problems do exploit, im- 
plicitly, the local behavior mentioned above. In particular, explicit schemes make use of 
polynomial approximations to the exponential whereas implicit schemes employ rational ap- 
proximations. It is well known that the stability properties of implicit schemes resulting 
from these approximations make their use more attractive than an explicit scheme. The 
price paid for this stability is the necessity to solve one or more large linear systems at each 
time step. On the other hand explicit schemes suffer from their requirement of using small 
time steps. The Krylov subspace approach presented in this paper avoids the solution ol 
large linear systems of equations necessitated by the use of implicit time integration schemes 
while retaining the stability properties of these schemes. The elimination of the necessity 
of solving large systems of equations makes the use of higher order polynomials and their 
concomitant higher order accuracy very attractive from the practical point of view. The 
immediate object of this work is to demonstrate the feasibility of using these methods to 
obtain steady state solutions of the Navier-Stokes equations. 

Krylov subspace methods form the basis of powerful linear system solvers such as the 
conjugate gradient and generalized minimum residual method. In section two of this paper 
we indicate how Krylov subspaces can also be used to approximate the product of a matrix 
exponential times a given vector. The basic idea has been suggested in several previous 
articles, see for example [3, 1, 9]. In the same section we also briefly discuss how the 
approximation can be used to solve time dependent nonlinear partial differential equations. 

Section 3 deals with the application of the method to the problem of incompressible vis- 
cous flow. The Navier-Stokes equations in primitive variables form are solved using a pseudo 
compressibility approach. In this work the spacial derivatives are approximated by second 
order centered differences on a staggered mesh. The resulting nonlinear differential equa- 
tions are then manipulated into a form suitable for the application of the Krylov exponential 
method. 
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2 The Krylov subspace approach 

In this section we address the problem of computing an approximation to a vector w = e~ A v 
by using polynomial approximation, i.e., we seek an approximation of the form: 

t~ A v « p m _i(A)v, (1) 

where p m - 1 is a polynomial of degree m— 1 . A systematic way in which such an approximation 
can be found is to start by observing that the vector w = e~ A v is an element of the Krylov 
subspace 

K m = span{v , Av, A m ~ 1 v}. 

One would like to get an element of this subspace that approximates e~ A v, ideally in some 
least squares sense. To facilitate working with vectors in K m , it is convenient to generate 
an orthonormal basis V m = [ui,i>2) v 3, • • • > v m ]. Starting with the vector uj = u/|]v ||2 we will 
generate the basis V m with the well-known Arnold! algorithm, described below. 


Algorithm: Arnoldi 

1. Initialize: 

Compute ui := u/||t/|| 2 . 

2. Iterate: Do j = 1,2, ...,m 

1. Compute w := Avj. 

2. Compute a set of j coefficients hij so that 

j 

w-w-Y, hijVi 
1=1 

is orthogonal to all previous v^s. 

3. Compute hj+tj := ||in ||2 and Vj + 1 := w/hj+ij. 

The above algorithm produces an orthonormal basis V m = [ui, v 2 , . . . , u m ], of the Krylov 
subspace K m . If we denote the m x m upper Hessenberg matrix consisting of the coefficients 
hij computed from the algorithm by II m we have the relation 

A Vfn — V m Hm “t (") 

Because of the orthogonality of the matrix V m we have the relation H m = V^AVm and as 
a result H m represents the projection of the linear transformation A to the subspace A m , 
with respect to the basis V m . The vector x opt = V m V^e~ A v is the projection of e~ A v on I< m , 
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i.e., it is the closest approximation to exp(— A)v from K m . Defining 0 = ||v|| 2 , h follows 
immediately that, 

V m V^e~ A v = 0V m Vlt~ A v x = V m Vle~ A V m t x . 

The optimal fit is therefore x op < = V m y op t in which y op t = 0V^e A V m e 1 . Unfortunately, 
y ov t is not practically computable, since it involves e~ A . The alternative which will be used 
throughout is to approximate V^e~ A V m by e~ Hm , leading to the approximation, y op t ~ 
0e~ Hm ei and, . 

e~ A v « 0V m e~ Hm t l . (3) 

We are now left with the problem of computing efficiently the vector e m ei which is 
similar to the problem we started with but typically of much smaller size. We will not discuss 
the aspect of this computation here. The reader is referred to [4, 3] for one approach based 
on rational approximation to the exponential. More standard techniques such as those based 
on spectral decomposition may also be used[2] although these are typically more expensive. 
More generally, one can approximate w = <f>(A)v for any function which is analytic around 
the eigenvalues of A in the same manner. 

The theory of the above approximation has been examined in Gallopoulos and Saad [3] 
and [9]. Applications to the solution of systems of Ordinary Differential Equations have been 
discussed by several authors [1, 5, 7, 8] 

We now give a brief overview on the ways in which the above approximation can be 
exploited to solve systems of Ordinary Differential Equations. We start with the simplest 
such equation, namely the homogeneous linear system 


w' = — Aw 


which can be solved by using a time integration scheme of the form 

w(t + At) = e~* tA w(t). 

We can exploit the previous approximation to evaluate the above expression at any given 
time t. Note that the accuracy of the above scheme is basically dictated by that of the 
approximation (3) and that if the solution is desired at a certain time then in theory one 
may be able to compute it in one single time-step provided m is sufficiently large to make 
the approximation accurate enough. 

To solve more realistic time dependent equations we must start by looking at the slightly 
more general case: 

w' = —Aw + r (4) 

in which we assume for now that r does not depend on time. The exact solution to the above 
equation is known to be 

w(t) = A _1 r + e~ tA (w 0 — A *r) (5) 

which can be rewritten as 

u>(f) = w 0 + t<f>(tA){r - Au> 0 )- (6) 
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where we define _ z 

m = 1 ^-. (") 

z 

Therefore, similarly to the homogeneous case consider above we can devise an integration 
scheme based on the following formula 

u;(t + At) = w(t) + At(j>(AtA)(r - Aw(t)). (8) 

The use of the above formula, as opposed to the more standard equation (5) avoids the 
solution of linear systems. The evaluation of <f>(H m )e i has been discussed in detail [9]. 

The simplest way to generalize the above formula to nonlinear equations of the form 

w + L(w ) = 0 (9) 

is to rewrite it in a given interval [t,t + A/] as 

w' = —Aw + (Aw — L(w )) = —Aw + r(u>(f)) 

in which A is the differential of L with respect to w at w(t ), and then to replace the variable 
term Aw — L(w) by some constant r c in the interval [i,t + At]. The original equation (9) is 
therefore replaced by an equation of the form (4) which is advanced from t to t + At by a 
formula of the form (8). Many variations to this basic scheme are possible. 

More generally, one can devise time marching procedures by attempting to exploit the 
expression for the exact solution of (9) at time t + At, which is given by 

w(t + At) = w(t) + [r(u>(f + t)) - Aw(t)] dr (10) 

Jo 

One can develop numerous schemes by applying quadrature formulas to the above equation, 
in the same spirit that prevails in the development of Runge Kutta methods. An important 
aspect of these is that the integrand involves the matrix exponential function. 

An important detail concerning the practical use of the above idea is that of estimating 
the error a posteriori. This is important because it allows us to determine whether or not the 
computation that has been performed using (3) is accurate enough. Such an estimator has 
been presented in [9]. In case the accuracy obtained is insufficient, one can easily backtrack 
by reducing the time step At. An important point here is that, according to the formula 
(3) the Krylov subspace need not be recomputed. We only need to redo the computation of 
exp(—AtH m )ei involving the small matrix. In addition it is known that roughly the error 
behaves like C x (At) m where C is a constant. This allows to determine an almost optimal 
At after the first attempt. In practice it is extremely rare that more than one backtracking 
step is performed which makes the process of adjusting time steps very economical. 

3 Application to Fluid Dynamics 

This section examines the use of the matrix exponential approach to solve computational 
fluid dynamics problems. In particular, the case of incompressible flow is addressed and 


5 


methods for steady state problems are presented. We begin with an introduction to the 
incompressible Navier-Stokes equations. The equations in nondimensional form are. 

K-^-A V + VP = -{V-V)V (11) 

Re 

V • V = 0 (12) 

In two dimensions, the velocity vector is made up of the x and y components of the velocity 

and P is the pressure. The gradient operator is V, V* is the divergence operator and A 
represents the Laplacian d^/dx^ + <9 2 /<9j/ 2 . The dissipation is scaled by Re, the Reynolds 
number. 


3.1 Steady State Case 

The artificial or pseudo compressibility method is the starting point in the application of 
exponential propagation by Krylov subspaces. The method involves adding a time derivative 
of pressure to the continuity equation and then solving the perturbed system. The equations 
to be solved are: 

V t -1-AV + (V- V)V = -VP (13) 

P T + /?V • V = 0 (14) 

where (3 is the compressibility parameter. For a mor detailed discussion of pseudo compress- 
ibility methods, the reader is refered to the early paper by Chorin [10] and the more recent 
work of Kwak and Chang [11, 12]. The usual method of solution of the above system is as 
follows: 

• knowing V n and P„ integrate equation (13) for V„+i 

• P n+1 =P„-Ar/?V-K + i 

Time integration of equation (13) can be achieved by either implicit or explicit methods. 


3.2 Use of Matrix Exponentials 


After discretization in space, the perturbed Navier-Stokes equations can be written as, 


+ f(W) = b. 


(15) 


In this system the vector W contains the unknown velocities and pressure at the nodes. The 
vector b contains the known boundary values, which are constant, and the function / which 
is nonlinear in W is defined by, 


m o = 


+ VP + (V • V)V 
£V • V 


(16) 
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The system (15) represents a nonlinear system of ordinary differential equations. 

In order to apply the methods of the previous section the nonlinear equations are lin- 
earized around the current value of W, W{t). 

W t + f(W(t)) + J(W - W(t)) + [f(W) - - J(W - W(*))] = b (17) 

After some manipulation we obtain the equation, 

W t + JW = b - f{W) + JW = r{W) (18) 

where J is the Jacobian of / evaluated at l V[t) or J = df /dW. 

The solution of equation (18) can be written in terms of matrix exponentials. 

W{t + At) = W(t) + e-< At - T > J [r(W(t + r)) - JW(t)) dr (19) 

Jo 

Application of this formula to advance the solution of equation (15) in time requires the 
approximation of the integral by a quadrature formula. This is complicated by the fact 
that the integrand is not an explicit function of r. However, the principle of approximately 
integrating equations similar to (19) using numerical quadrature has been at the basis of 
methods of Runge Kutta type. 


3.3 Solution of the Driven Cavity Problem 

The problem of viscous flow in a lid driven cavity was used as a test case for solutions 
of the Navier-Stokes equations via equation (19). This problem has the attractive feature 
of no inflow or outflow boundaries. The boundary conditions for this case are simply the 
specification of the velocities on the boundaries. Since the problem is discretized on a 
staggered mesh, no pressure boundary condition is required. 

The integral in equation (19) was evaluated in two ways. The evaluation of r at W(t) 
yields, after evaluation of the integral, the first of two formula for evaluating W(t + At). 

W(t + At) = W{t) + A t(I - e- AU )(AtJ)-\b - f(W(t))) (20) 


We call this formula the left point rule because the value of the term in brackets in equa- 
tion (19) is assumed to have a constant value in the entire interval from t to t + At. This 
value is the value of the term at the left of the interval. The midpoint rule is arrived at by 
similar means. Evaluating the term in brackets at the midpoint of the interval results in the 
following formula for W{t - f At). 


W(t + At) = W(t) + Ate~^ J 


r(W{t + ^)) - JW(t) 


( 21 ) 


The left point rule is used to predict the value of r(W(t + ^)) and then the mid point 
formula is used to correct the solution. 

Both the left point method alone and the predictor-corrector method composed of the 
left point and mid point rules were validated on the driven cavity problem. An equally 
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spaced staggered grid was used for space discretization. The methods of the previous section 
were used to evaluate the matrix exponential terms which appear in the two formulas. The 
matrix exponential methods shown above required only the evaluation of a matrix vector 
product. In this particular case finite differences were used to approximate the Jacobian 
matrix vector product. 

+ en) - 


Jv = 


( 22 ) 


This technique makes the explicit computation of the matrix J and its storage unnecessary. 

Converged solutions for both the methods were obtained on meshes with 32 points in 
each of the coordinate directions. The Reynolds numbers for both tests was 400 and the 
compressibility parameter /? and time step were 2.5 and 0.01 respectively. The initial con- 
ditions were velocity and pressure set to zero inside the cavity. A plot of the vorticity of 
the converged solution is shown in figure 2. Both methods gave identical solutions. Only 
the results of the mid point formula are shown. Figure 1 contains three plots of convergence 
history for the first 2000 iterations of the method. The three plots represent the behavior of 
the logrithm of the norm of the vector b- f(W). The information is separated into the parts 
of the vector corresponding to the steady state u momentum, v momentum and continuity 
equations. The continuity equation data is a direct measure of the divergence of the velocity 
field. All the plots show convergence of the algorithm. 
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Figure 2: Vorticity 
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4 Conclusion 


The results shown in this paper, although preliminary, show the feasibility of using matrix 
exponential methods based on Krylov subspaces, for solving fluid dynamics problems. The 
schemes that we used for handling the integral are based on the simplest quadrature rules 
possible. Yet the results obtained show good accuracy and the corresponding schemes are 
fairly efficient, ft has been encouraging to see that the method seems very insensitive to 
the choice of the time step suggesting robust stability properties. A full study of the exact 
stability behavior remains to be done. We note that a number of related schemes presented 
in [3] for linear problems and using a different method for handling the integral, have been 
examined from this point of view and were shown to be either stable or Aq stable. 

It is hoped that more sophisticated integration schemes will lead to efficient procedures 
similar to Runge Kutta schemes but based on exploiting the matrix exponential. 
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